Load packages and configure the grid resolution and source filename:
library(dplyr)
library(ggplot2)
library(dggridR)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)
library(prettygrids)
library(geosphere)
library(stringr)
world <- ne_countries(scale = "medium", returnclass = "sf")
resolution <- 9
filename <- "occurrence_minimal_20200212"
Use the dggridR package to create a discrete global grid, and fetch the grid cell surface area to be included in the map title:
dggs <- dgconstruct(projection = "ISEA", res = resolution, resround = "up")
area <- dggetres(dggs)$area_km[resolution + 1]
area_str <- format(signif(area, 2), big.mark = ",")
deg_per_km <- destPoint(c(0, 0), 90, 1000)[1,1]
square_side_km <- sqrt(area)
square_side_deg <- square_side_km * deg_per_km
square_side_deg_str <- format(signif(square_side_deg, 2), big.mark = ",")
This notebook uses a A CSV file with occurrence records exported from the database database. Alternatively the robis R package can be used, but downloading the entire database can take a while. This downloads the zipped CSV file and aggregates the data into locations with number of records:
if (!file.exists("grouped.temp")) {
if (!file.exists("data.zip")) {
download.file(paste0("https://download.obis.org/export/", filename, ".zip"), "data.zip")
}
df_raw <- read.csv(unz("data.zip", paste0(filename, ".csv")), stringsAsFactors = FALSE)
df_grouped <- df_raw %>%
group_by(decimallongitude = round(decimallongitude, 2), decimallatitude = round(decimallatitude, 2), speciesid) %>%
summarize(records = n(), species = length(unique(speciesid)))
save(df_grouped, file = "grouped.temp")
file.remove("data.zip")
} else {
load("grouped.temp")
}
Next we assign cell IDs to each location and caluclate statistics per cell:
df_grouped$cell <- dgtransform(dggs, df_grouped$decimallatitude, df_grouped$decimallongitude)
stats <- df_grouped %>%
group_by(cell) %>%
summarise(records = sum(records), species = length(unique(speciesid)))
map_theme <- theme(
axis.ticks.length = unit(0, "pt"),
line = element_blank(), rect = element_blank(),
axis.text = element_blank(), axis.title = element_blank()
)
map_scale_records <- scale_fill_viridis_c(
option = "plasma",
trans = "log10",
name = "Records",
breaks = c(1, 10, 100, 1000, 10000),
labels = function(x) format(x, scientific = FALSE),
limits = c(1, 100000),
oob = scales::squish
)
map_scale_species <- scale_fill_viridis_c(
option = "plasma",
trans = "log10",
name = "Species",
labels = function(x) format(x, scientific = FALSE)
)
make_subtitle <- function(type, projection, marineregions = FALSE) {
s <- paste0("Number of ", type, " per ISEA3H grid cell of approximately ", area_str, " square km (this corresponds to a ", square_side_deg_str, " by ", square_side_deg_str, " degree cell at the equator).\n", projection, ". Ocean Biogeographic Information System (OBIS), 2020.")
if (marineregions) {
s <- paste0(s, " Maritime boundaries from marineregions.org.")
}
return(s)
}
Depending on the projection the grids produced by dggridR can cause artefacts in your map. Here we are using the prettygrids::make_grid utility function (still under development) to fix some of these problems.
grid_pac <- make_grid(offset = -180, res = resolution) %>%
merge(stats, by.x = "cell", by.y = "cell")
crs <- "+proj=robin +lon_0=180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
bbox <- st_bbox(st_as_sfc(st_bbox(c(xmin = 18, xmax = 248, ymin = -55, ymax = 55), crs = 4326)) %>% st_transform(crs))
ggplot() +
geom_sf(data = grid_pac %>% st_transform(crs), aes(fill = records), size = 0, color = NA) +
geom_sf(data = world %>% prepare_grid(offset = -180) %>% st_transform(crs), fill = "#cccccc", color = "#999999", size = 0.1) +
coord_sf(crs = crs, xlim = c(bbox$xmin, bbox$xmax), ylim = c(bbox$ymin, bbox$ymax)) +
map_theme +
map_scale_records +
labs(
title = "Number of records",
subtitle = make_subtitle("records", "Pacific centered Robinson projection")
)
ggplot() +
geom_sf(data = grid_pac %>% st_transform(crs), aes(fill = species), size = 0, color = NA) +
geom_sf(data = world %>% prepare_grid(offset = -180) %>% st_transform(crs), fill = "#cccccc", color = "#999999", size = 0.1) +
coord_sf(crs = crs, xlim = c(bbox$xmin, bbox$xmax), ylim = c(bbox$ymin, bbox$ymax)) +
map_theme +
map_scale_species +
labs(
title = "Number of species",
subtitle = make_subtitle("species", "Pacific centered Robinson projection")
)
bbox_usa <- st_bbox(st_as_sfc(st_bbox(c(xmin = 130, xmax = 297, ymin = -30, ymax = 70), crs = 4326)) %>% st_transform(crs))
usa <- st_read("eez_v11_0_360_dissolved_simplified005.gpkg", crs = 4326) %>%
filter(SOVEREIGN1 == "United States")
## Reading layer `eez_v11_0_360_dissolved_simplified005' from data source `/Users/pieter/Desktop/dev/notebooks/density_maps/eez_v11_0_360_dissolved_simplified005.gpkg' using driver `GPKG'
## Simple feature collection with 158 features and 31 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 0 ymin: -85.5625 xmax: 360 ymax: 86.98832
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
ggplot() +
geom_sf(data = grid_pac %>% st_transform(crs), aes(fill = records), size = 0, color = NA) +
geom_sf(data = world %>% prepare_grid(offset = -180) %>% st_transform(crs), fill = "#cccccc", color = "#999999", size = 0.1) +
geom_sf(data = usa %>% st_transform(crs), fill = NA, color = "#000000", size = 1) +
coord_sf(crs = crs, xlim = c(bbox_usa$xmin, bbox_usa$xmax), ylim = c(bbox_usa$ymin, bbox_usa$ymax)) +
map_theme +
map_scale_records +
labs(
title = "Number of records",
subtitle = make_subtitle("records", "Pacific centered Robinson projection", marineregions = TRUE)
)
ggplot() +
geom_sf(data = grid_pac %>% st_transform(crs), aes(fill = species), size = 0, color = NA) +
geom_sf(data = world %>% prepare_grid(offset = -180) %>% st_transform(crs), fill = "#cccccc", color = "#999999", size = 0.1) +
geom_sf(data = usa %>% st_transform(crs), fill = NA, color = "#000000", size = 1) +
coord_sf(crs = crs, xlim = c(bbox_usa$xmin, bbox_usa$xmax), ylim = c(bbox_usa$ymin, bbox_usa$ymax)) +
map_theme +
map_scale_species +
labs(
title = "Number of species",
subtitle = make_subtitle("species", "Pacific centered Robinson projection", marineregions = TRUE)
)
north <- st_set_crs(st_as_sf(as(raster::extent(-180, 180, 0, 90), "SpatialPolygons")), 4326)
grid_ortho <- make_grid(offset = NA, res = resolution) %>%
merge(stats, by.x = "cell", by.y = "cell") %>%
filter(st_intersects(geometry, north, sparse = FALSE))
crs <- "+proj=laea +lat_0=52 +lon_0=0 +x_0=4321000 +y_0=3210000 +datum=WGS84 +units=m +no_defs"
bbox <- st_sfc(st_polygon(list(matrix(c(-140, 40, -60, 40, 60, 40, 140, 40, -140, 40), ncol = 2, byrow = TRUE))), crs = 4326)
bbox_trans <- bbox %>% st_transform(crs)
ggplot() +
geom_sf(data = grid_ortho %>% st_transform(crs), aes(fill = records), size = 0, color = NA) +
geom_sf(data = world %>% st_transform(crs), fill = "#cccccc", color = "#999999", size = 0.1) +
coord_sf(crs = crs,
xlim = c(min(st_coordinates(bbox_trans)[,1]), max(st_coordinates(bbox_trans)[,1])),
ylim = c(min(st_coordinates(bbox_trans)[,2]), max(st_coordinates(bbox_trans)[,2]))) +
map_theme +
map_scale_records +
labs(
title = "Number of records",
subtitle = make_subtitle("records", "Lambert azimuthal equal-area projection")
)
ggplot() +
geom_sf(data = grid_ortho %>% st_transform(crs), aes(fill = species), size = 0, color = NA) +
geom_sf(data = world %>% st_transform(crs), fill = "#cccccc", color = "#999999", size = 0.1) +
coord_sf(crs = crs,
xlim = c(min(st_coordinates(bbox_trans)[,1]), max(st_coordinates(bbox_trans)[,1])),
ylim = c(min(st_coordinates(bbox_trans)[,2]), max(st_coordinates(bbox_trans)[,2]))) +
map_theme +
map_scale_species +
labs(
title = "Number of species",
subtitle = make_subtitle("species", "Lambert azimuthal equal-area projection")
)